library(cregg)
library(ggplot2)
library(texreg)


d <- read.csv("NI_border_replication_data.csv", stringsAsFactors=T)

#### data management ####
# relevel
d$hrprotection <- relevel(d$hrprotection, ref =  "New international tribunal" ) 
#d$hrprotection <- relevel(d$hrprotection, ref =  "UK only" )
d$public_spending <- relevel(d$public_spending, ref = "Same")
d$location <- relevel(d$location, ref =  "North/South")
d$executive <- relevel(d$executive, ref = "All parties")
d$legislation <- relevel(d$legislation, ref =  "Forced coalition")


# set feature labels
attr(d$hrprotection, "label") <- "COURTS"
attr(d$public_spending, "label") <- "ECONOMIC COMPENSATION"
attr(d$location, "label") <- "TERRITORIAL ISSUES"
attr(d$executive, "label") <- "EXECUTIVE"
attr(d$legislation, "label") <- "LEGISLATION"
attr(d$group3, "label") <- "Community"

# set feature labels for controls
d$income_rec <- factor(d$income_rec, levels = c('?0 to ?14,999 per year', '?15,000 to ?34,999 per year', '?35,000 to ?59,999 per year', '?60,000 to ?99,999 per year', '?100,000+ per year'))
table(d$income_rec)
d$education_lab <- factor(d$education_lab, levels = c("No formal qualifications", 'CSE grade 1, GCE O level, GCSE, School + CSE grades 2-5', 'GCE A level or Higher Certificate + Irish Leaving certificate', 'City & Guilds certificates + other qualifications', 'University or CNAA first degree', 'University or CNAA higher degree', 'PhD'))
attr(d$age, "label") <- "AGE"
attr(d$sex, "label") <- "SEX"
attr(d$education_lab, "label") <- "EDUCATION"
attr(d$income_rec, "label") <- "INCOME"

#### AMCE ####
# without weights

m <- selected ~ location + executive + legislation  + public_spending +  hrprotection

amces <- cj(d, m, id = ~id,  by = ~ group3)
head(amces[c("feature", "level", "estimate", "std.error")], 10L) 

# plot 
plot(amces, group = 'group3',  header_fmt = "%s") +
  scale_color_manual(values = c("#339900", "#FFCC00", "#FF6633"), na.translate = F) +
  facet_wrap(~feature, ncol = 1L,
             scales = "free_y", strip.position = "right") +
  theme(legend.title = element_blank(),
        strip.text.y = element_blank()) 


#### MARGINAL MEANS ####
marginalmeans <- cj(d, m, id = ~id, by = ~ group3, estimate = 'mm')
head(marginalmeans[c("feature", "level", "estimate", "std.error")], 20L) 

# plot 
plot(marginalmeans, group = 'group3', vline = 0.5,  header_fmt = "%s") +
  scale_color_manual(values = c("#339900", "#FFCC00", "#FF6633"), na.translate = F) +
  scale_x_continuous(limits = c(.1, .9), breaks = c(.1, .3, .5, .7, .9)) +
  facet_wrap(~feature, ncol = 1L,
             scales = "free_y", strip.position = "right") +
  theme(legend.title = element_blank(),
        strip.text.y = element_blank()) 

# APPENDIX
#### Regression table for AMCE ####
# nationalists
m_amce_nat <- lm(selected ~ location + executive + legislation  + public_spending +  hrprotection, subset(d, group3=='Nationalists'))
summary(m_amce_nat)

# unionists
m_amce_uni <- lm(selected ~ location + executive + legislation + public_spending + hrprotection, subset(d, group3=='Unionists'))
summary(m_amce_uni)

# neither
m_amce_neither <- lm(selected ~ location + executive + legislation + public_spending + hrprotection, subset(d, group3=='Neither'))
summary(m_amce_neither)

# create table 
wordreg(list(m_amce_nat, m_amce_uni, m_amce_neither),
        file = "NI_amce.docx", ci.force = TRUE,
        groups = list("TERRITORIAL ISSUES" = 2, "EXECUTIVE" = 3:5, "LEGISLATION" = 6:7, "ECONOMIC COMPENSATION" = 8:9, "COURTS" = 10:12),
        custom.coef.names = c("(Intercept)", "East/West", "1/4 from each", "Largest Unionist+Nationalist", "Majority rule", "1/4 from both", "Majority of members", "10%", "5%", "European and international courts", "New international tribunal", "UK, EU and international courts"),
        custom.model.names = c("Nationalists","Unionists", "Neither"))
unlink("NI_amce.doc")

#### Regression table for MM ####
mm <- marginalmeans[c("feature", "level", "estimate", "p", "lower", "upper", "group3")]
# scientific is an issue!
write.table(format(mm, digits=2), file = "NI_mm.txt", sep = ";", quote = FALSE, row.names = F)

#### AMCE with weights ####
# with weights, by community 3 ways
amces_weights <- cj(d, m, id = ~id,  weights = ~weight, by = ~ group3)
head(amces_weights[c("feature", "level", "estimate", "std.error")], 10L) 

# plot 
plot(amces_weights, group = 'group3', header_fmt = "%s") +
  scale_color_manual(values = c("#339900", "#FFCC00", "#FF6633"), na.translate = F) +
  # xlim(-.4, .4) +
  facet_wrap(~feature, ncol = 1L,
             scales = "free_y", strip.position = "right") +
  theme(legend.title = element_blank(),
        strip.text.y = element_blank()) 

#### MM with weights ####
marginalmeans_weights <- cj(d, m, id = ~id,  weights = ~weight, by = ~ group3, estimate = 'mm')
head(marginalmeans_weights[c("feature", "level", "estimate", "std.error")], 20L) 

# plot 
plot(marginalmeans_weights, group = 'group3', vline = 0.5,  header_fmt = "%s") +
  scale_color_manual(values = c("#339900", "#FFCC00", "#FF6633"), na.translate = F) +
  scale_x_continuous(limits = c(.1, .9), breaks = c(.1, .3, .5, .7, .9)) +
  facet_wrap(~feature, ncol = 1L,
             scales = "free_y", strip.position = "right") +
  theme(legend.title = element_blank(),
        strip.text.y = element_blank()) 

#### AMCE with controls ####
# without weights
# formula model conjoint with controls
m_controls <- selected ~ location + executive + legislation  + public_spending +  hrprotection  + age + sex + income_rec + education_lab

amces_controls <- cj(d, m_controls, id = ~id, by = ~ group3)

class(d$hrprotection)
table(d$income_rec)

# plot 
plot(amces_controls, group = 'group3', header_fmt = "%s") +
  scale_color_manual(values = c("#339900", "#FFCC00", "#FF6633"), na.translate = F) +
  xlim(-.4, .4) +
  facet_wrap(~feature, ncol = 1L,
             scales = "free_y", strip.position = "right") +
  theme(legend.title = element_blank(),
        strip.text.y = element_blank()) 

#### Regression table for AMCE with controls without weights ####
# nationalists
m_amce_contr_nat <- lm(selected ~ location + executive + legislation + public_spending + hrprotection  + age + sex + income_rec + education_lab, subset(d, group3=='Nationalists'))
summary(m_amce_contr_nat)

# unionists
m_amce_contr_uni <- lm(selected ~ location + executive + legislation + public_spending + hrprotection  + age + sex + income_rec + education_lab, subset(d, group3=='Unionists'))
summary(m_amce_contr_uni)

# neither
m_amce_contr_neither <- lm(selected ~ location + executive + legislation + public_spending + hrprotection  + age + sex + income_rec + education_lab, subset(d, group3=='Neither'))
summary(m_amce_contr_neither)

# create table 
wordreg(list(m_amce_contr_nat, m_amce_contr_uni, m_amce_contr_neither), 
        file = "NI_amce_controls.doc", ci.force = TRUE, 
        groups = list("TERRITORIAL ISSUES" = 2, "EXECUTIVE" = 3:5, "LEGISLATIVE" = 6:7, "ECONOMIC COMPENSATION" = 8:9, "COURTS" = 10:12, "AGE" = 13:15, "SEX" = 16, "INCOME" = 17:20, "EDUCATION" = 21:26),
        custom.coef.names = c("(Intercept)", "East/West", "1/4 from each", "Largest Unionist+Nationalist", "Majority rule", "1/4 from both", "Majority of members", "10%", "5%", "European and international courts", "New international tribunal", "UK, EU and international courts", "AGE: 25-44 years", "45-64 years", "65+ years", "Male", "?15,000 to ?34,999 per year", "?35,000 to ?59,999 per year", "?60,000 to ?99,999 per year", "?100,000+ per year", "CSE grade 1, GCE O level, GCSE, School + CSE grades 2-5", "GCE A level or Higher Certificate + Irish Leaving certificate", "City & Guilds certificates + other qual.", "University or CNAA first degree", "University or CNAA higher degree", "PhD"),
        custom.model.names = c("Nationalists","Unionists", "Neither"))
unlink("NI_amce_controls.doc")

#### AMCE first conjoint only ####
d_first <- d[d$task==1,]
# set feature labels
attr(d_first$hrprotection, "label") <- "COURTS"
attr(d_first$public_spending, "label") <- "ECONOMIC COMPENSATION"
attr(d_first$location, "label") <- "TERRITORIAL ISSUES"
attr(d_first$executive, "label") <- "EXECUTIVE"
attr(d_first$legislation, "label") <- "LEGISLATIVE"
attr(d_first$group3, "label") <- "Community"

amces_first <- cj(d_first, m, id = ~id,  weights = ~weight, by = ~ group3, header_fmt = "%s")
plot(amces_first, group = 'group3') +
  scale_color_manual(values = c("#339900", "#FFCC00", "#FF6633"), na.translate = F) +
  xlim(-.6, .6) + # LIMITS NEED CHANGING?
  facet_wrap(~feature, ncol = 1L,
             scales = "free_y", strip.position = "right") +
  theme(legend.title = element_blank(),
        strip.text.y = element_blank()) 

rm(d_first, amces_first)

#### AMCE last conjoint only ####
d_last <- d[d$task==5,]
# set feature labels
attr(d_last$hrprotection, "label") <- "COURTS"
attr(d_last$public_spending, "label") <- "ECONOMIC COMPENSATION"
attr(d_last$location, "label") <- "TERRITORIAL ISSUES"
attr(d_last$executive, "label") <- "EXECUTIVE"
attr(d_last$legislation, "label") <- "LEGISLATIVE"
attr(d_last$group3, "label") <- "Community"

amces_last <- cj(d_last, m, id = ~id,  weights = ~weight, by = ~ group3)
plot(amces_last, group = 'group3', header_fmt = "%s") +
  scale_color_manual(values = c("#339900", "#FFCC00", "#FF6633"), na.translate = F) +
  xlim(-.6, .6) +
  facet_wrap(~feature, ncol = 1L,
             scales = "free_y", strip.position = "right") +
  theme(legend.title = element_blank(),
        strip.text.y = element_blank()) 

#### Omnibus F test ####
# cj_anova does not support weights
cj_diff <- cj_anova(d, m, id = ~id, by = ~ group3)
cj_diff 
write.table(format(cj_diff, digits=2), file = "NI_omnibus_f_test.txt", sep = ",", quote = FALSE, row.names = F)






